/*
 * gsl_app.cpp
 *
 *  Created on: May 27, 2011
 *      Author: nzhao
 */

#include "share/Usrdef.h"
#include "share/gsl_app.h"


ostream             &operator << (ostream &c, gsl_complex z)
                    {
                        c << "(" << GSL_REAL(z) << ", ";
                        c        << GSL_IMAG(z) << ")";
                        return c;
                    };

gsl_complex         operator - (gsl_complex z)
                    {return gsl_complex_negative(z);};

gsl_complex         operator + (gsl_complex z1, gsl_complex z2)
                    {return gsl_complex_add(z1,z2);};
gsl_complex         operator - (gsl_complex z1, gsl_complex z2)
                    {return gsl_complex_sub(z1,z2);};
gsl_complex         operator * (gsl_complex z1, gsl_complex z2)
                    {return gsl_complex_mul(z1,z2);};
gsl_complex         operator / (gsl_complex z1, gsl_complex z2)
                    {return gsl_complex_div(z1,z2);};

gsl_complex         operator + (double d, gsl_complex z)
                    {return gsl_complex_add_real(z, d);};
gsl_complex         operator + (gsl_complex z, double d)
                    {return gsl_complex_add_real(z, d);};

gsl_complex         operator - (double d, gsl_complex z)
                    {return gsl_complex_add_real(-z, d);};
gsl_complex         operator - (gsl_complex z, double d)
                    {return gsl_complex_sub_real(z, d);};

gsl_complex         operator * (double d, gsl_complex z)
                    {return gsl_complex_mul_real(z,d);};
gsl_complex         operator * (gsl_complex z, double d)
                    {return gsl_complex_mul_real(z,d);};

gsl_complex         operator / (double d, gsl_complex z)
                    {return d*gsl_complex_inverse(z);};
gsl_complex         operator / (gsl_complex z, double d)
                    {return gsl_complex_div_real(z,d);};

gsl_complex gsl_double_2_complex (double d)
{
	gsl_complex res;
	GSL_SET_COMPLEX(&res, d, 0.);
	return res;
}

int gsl_vector_double_2_complex (gsl_vector * v, gsl_vector_complex * vc)
{
	int i;
	gsl_complex z;
	int n;

	n = (int)(v->size);
//	vc = gsl_vector_complex_alloc(n);
	for (i = 0; i < n; i++)
	{
		z = gsl_double_2_complex(gsl_vector_get(v, i));
		gsl_vector_complex_set(vc, i, z);
	}
	return 1;
}

int gsl_matrix_double_2_complex (gsl_matrix * m, gsl_matrix_complex * mc)
{
	int i, j;
	gsl_complex z;
	int n1, n2;

	n1 = (int)(m->size1);
	n2 = (int)(m->size2);

//	mc = gsl_matrix_complex_alloc(n1,n2);
	for (i = 0; i < n1; i++)
	{
	    for(j = 0; j < n2; j++)
	    {
            z = gsl_double_2_complex(gsl_matrix_get(m, i, j));
            gsl_matrix_complex_set(mc, i, j, z);
	    }
	}
	return 1;
}

int gsl_matrix_complex_conj_transpose(gsl_matrix_complex * m)
{
	int i, j, n1, n2;
	n1 = (int)(m->size1);
	n2 = (int)(m->size2);
	gsl_complex z;

	gsl_matrix_complex_transpose(m);
	for (i = 0; i < n1; i++)
	{
		for(j = 0; j < n2; j++)
		{
			z = gsl_matrix_complex_get(m, i, j);
			gsl_matrix_complex_set(m, i, j, gsl_complex_conjugate(z));
		}
	}
	return 1;
}

double gsl_matrix_trace(gsl_matrix * m)
{
	double res;
	size_t dim, i;

	dim = m->size1;


	res = 0.;
	for (i = 0; i < dim; i++)
	{
		res +=  gsl_matrix_get(m, i, i);
	}
	return res;
}


gsl_complex gsl_matrix_complex_trace(gsl_matrix_complex * m)
{
	gsl_complex res;
	size_t dim, i;

	dim = m->size1;


	GSL_SET_COMPLEX (&res, 0., 0.);
	for (i = 0; i < dim; i++)
	{
		res = res + gsl_matrix_complex_get(m, i, i);
	}
	return res;
}
gsl_complex gsl_matrix_complex_trace1(gsl_matrix_complex * m, double * coeff)
{
	gsl_complex res;
	size_t dim, i;

	dim = m->size1;


	GSL_SET_COMPLEX (&res, 0., 0.);
	for (i = 0; i < dim; i++)
	{
		res = res + (gsl_matrix_complex_get(m, i, i) * coeff[i]);
	}
	return res;
}

int gsl_matrix_row_transform(size_t nrow, double * val, gsl_matrix * m)
{
	size_t n1, n2, i, j;

	n1 = m->size1;
	n2 = m->size2;

	if (n1!=nrow)
	{
		GSL_ERROR("dim not match", GSL_ERANGE);
	}

	double v_data_j, vi;
	gsl_vector_view v_row;
	for (i = 0; i < nrow; i++)
	{
		v_row = gsl_matrix_row(m, i);
		vi = val[i];
		for (j = 0; j < n2; j++)
		{
			v_data_j = gsl_vector_get(&v_row.vector, j);
			gsl_vector_set(&v_row.vector, j, v_data_j * vi);
		}
	}
	return 1;
}

int gsl_matrix_complex_row_transform(size_t nrow, gsl_complex * val, gsl_matrix_complex * m)
{
	size_t n1, n2, i, j;

	n1 = m->size1;
	n2 = m->size2;

	if (n1!=nrow)
	{
		GSL_ERROR("dim not match", GSL_ERANGE);
	}

	gsl_complex v_data_j, vi;
	gsl_vector_complex_view v_row;
	for (i = 0; i < nrow; i++)
	{
		v_row = gsl_matrix_complex_row(m, i);
		vi = val[i];
		for (j = 0; j < n2; j++)
		{
			v_data_j = gsl_vector_complex_get(&v_row.vector, j);
			gsl_vector_complex_set(&v_row.vector, j, v_data_j * vi);
		}
	}
	return 1;
}
